% This function finds the wage (w), the weight of imported final good in preferences (varrho)
% and the transfers measured in domestic labor (Timplied, which correspond to T/w) that 
% generate a desired expenditure share in non-tradables (param.shareN) and final
% good import share (param.shareM), given firm behavior (sD,sX).

% For consistency, it also computes the expenditure share in the imported final good (sM)
% and the expenditure share in the non-traded good (sN).

% See Calibration Appendix

function [w,sM,sN, varrho,Timplied] = calibrateaNM_byHand(sD, sX, phi,fc,fx,param)

eta = param.eta;
g = param.gamma;
eps = param.eps;
sg = param.sg;
beta = param.beta;
qD = param.qD;
theta = param.theta;
vartheta = param.vartheta;
zeta = param.zeta;
pM = param.pM;
AN = param.AN;
FXM = param.FXM; FX = param.FX; FM = param.FM;

% Constants

k1 =  qD^(g*(sg-1)) * (1-g)^((1-g)*(sg-1)) * g^(g*(sg-1)) * g * (sg-1)^sg/sg^sg * eta; % kappa2_wiggle
k2 = qD^(g*(sg-1)) * (sg-1)^(sg-1)/sg^sg * (1-g)^((1-g)*(sg-1))* g^(g*(sg-1)); % kappa3
k3 = k2 * (theta-1)/theta; % kappa3*times (theta-1)/theta
k5 = (sg / (sg-1))^(1-sg) * (1-g)^((1-g)*(sg-1)) * (g*beta^(eps/(eps-1))*qD)^(g*(sg-1)); % kappa1
k6 = mean(k1*(sD.^(-1)-1).^(1/(eps-1)/eta) .* fc + k2 .* ((sD<1).*(sX==0)) * FM + k3 * (sX./(1-sX)).^(theta/(theta-1)) .*fx + k2 * ((sX>0).*(sD==1)) * FX + k2 * ((sX>0).*(sD<1)) * FXM ); % \varPsi
k7 = mean(phi.^(sg-1) .* (sX./(1-sX)) .*sD.^(-g*(sg-1)/(eps-1))); % \Upsilon
k8 = mean(phi.^(sg-1) .* (1./(1-sX)) .*sD.^(1-g*(sg-1)/(eps-1))); %  \Gamma

P = (sg/ (sg-1))^(1/(1-g)) * (1-g)^(-1) * (g*beta^(eps/(eps-1))*qD)^(-g/(1-g)) * (mean(phi.^(sg-1).*sD.^(-g*(sg-1)/(eps-1))))^(1/((1-sg)*(1-g))); % This is P/w
bundle = (1-zeta)/zeta / (1-vartheta);
varrho  = 1 - (1/param.shareN-1)^bundle * (param.aN/(1-param.aN))^(vartheta*bundle) * (AN*P)^((zeta-1)/zeta) * (1- param.shareM)^(1/zeta);
w = pM / P * ( ((1-varrho)/varrho)^zeta * param.shareM/(1- param.shareM))^(1/(zeta-1));
sM = 1 / (1 + ((1-varrho)/varrho)^zeta * (P/(pM/w))^(1-zeta));
sN = 1 / (1 + ((1-param.aN)/param.aN)^vartheta * (AN * (varrho^zeta * (pM/w)^(1-zeta) + (1-varrho)^zeta * P^(1-zeta))^(1/(1-zeta))  )^(1-vartheta) );

Timplied = ( 1/(1-sM) * ( 1 - (sg-1)/sg * g *k5 *k8 * P^((sg-1)*(1-g)) ) - (1-(sg-1)/sg*g) * (  k5 * k7 * P^((sg-1)*(1-g)) + 1) ) / ( (sg-1)/sg*(1-g)*(k5 * k7 * P^((sg-1)*(1-g))+1) +  k6*P^((sg-1)*(1-g)) + sN/ ((1-sM)*(1-sN)) * (1- g* (sg-1)/sg * k5*k8* P^((sg-1)*(1-g))) ) * param.L; 
end